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Disclaimer 


The  US  Department  of  Commerce  makes  no  warranty,  expressed  or  implied,  to  users  of  the  Fire  Dy- 
namics Simulator  (FDS),  and  accepts  no  responsibility  for  its  use.  Users  of  FDS  assume  sole  responsibility 
under  Federal  law  for  determining  the  appropriateness  of  its  use  in  any  particular  application;  for  any  con- 
clusions drawn  from  the  results  of  its  use;  and  for  any  actions  taken  or  not  taken  as  a result  of  analyses 
performed  using  these  tools. 

Users  are  warned  that  FDS  is  intended  for  use  only  by  those  competent  in  the  fields  of  fluid  dynamics, 
thermodynamics,  combustion,  and  heat  transfer,  and  is  intended  only  to  supplement  the  informed  judgment 
of  the  qualified  user.  The  software  package  is  a computer  model  that  may  or  may  not  have  predictive 
capability  when  applied  to  a specific  set  of  factual  circumstances.  Lack  of  accurate  predictions  by  the  model 
could  lead  to  erroneous  conclusions  with  regard  to  fire  safety.  All  results  should  be  evaluated  by  an  informed 
user. 

Throughout  this  document,  the  mention  of  computer  hardware  or  commercial  software  does  not  con- 
stitute endorsement  by  NIST,  nor  does  it  indicate  that  the  products  are  necessarily  those  best  suited  for  the 
intended  purpose. 
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1 Introduction 


The  idea  that  the  dynamics  of  a fire  might  be  studied  numerically  dates  back  to  the  beginning  of  the  com- 
puter age.  Indeed,  the  fundamental  conservation  equations  governing  fluid  dynamics,  heat  transfer,  and 
combustion  were  first  written  down  over  a century  ago.  Despite  this,  practical  mathematical  models  of  fire 
(as  distinct  from  controlled  combustion)  are  relatively  recent  due  to  the  inherent  complexity  of  the  problem. 
Indeed,  in  his  brief  history  of  the  early  days  of  fire  research.  Hoyt  Hottel  noted  “A  case  can  be  made  for  fire 
being,  next  to  the  life  processes,  the  most  complex  of  phenomena  to  understand"  [1], 

The  difficulties  revolve  about  three  issues:  First,  there  are  an  enormous  number  of  possible  fire  scenarios 
to  consider  due  to  their  accidental  nature.  Second,  the  physical  insight  and  computing  power  necessary  to 
perform  all  the  necessary  calculations  for  most  fire  scenarios  are  limited.  Any  fundamentally  based  study 
of  fires  must  consider  at  least  some  aspects  of  bluff  body  aerodynamics,  multi-phase  flow,  turbulent  mixing 
and  combustion,  radiative  transport,  and  conjugate  heat  transfer;  all  of  which  are  active  research  areas  in 
their  own  right.  Finally,  the  “fuel”  in  most  fires  was  never  intended  as  such.  Thus,  the  mathematical  models 
and  the  data  needed  to  characterize  the  degradation  of  the  condensed  phase  materials  that  supply  the  fuel 
may  not  be  available.  Indeed,  the  mathematical  modeling  of  the  physical  and  chemical  transformations  of 
real  materials  as  they  bum  is  still  in  its  infancy. 

In  order  to  make  progress,  the  questions  that  are  asked  have  to  be  greatly  simplified.  To  begin  with, 
instead  of  seeking  a methodology  that  can  be  applied  to  all  fire  problems,  we  begin  by  looking  at  a few 
scenarios  that  seem  to  be  most  amenable  to  analysis.  Hopefully,  the  methods  developed  to  study  these  “sim- 
ple” problems  can  be  generalized  over  time  so  that  more  complex  scenarios  can  be  analyzed.  Second,  we 
must  learn  to  live  with  idealized  descriptions  of  fires  and  approximate  solutions  to  our  idealized  equations. 
Finally,  the  methods  should  be  capable  of  systematic  improvement.  As  our  physical  insight  and  computing 
power  grow  more  powerful,  the  methods  of  analysis  can  grow  with  them. 

To  date,  three  distinct  approaches  to  the  simulation  of  fires  have  emerged.  Each  of  these  treats  the 
fire  as  an  inherently  three  dimensional  process  evolving  in  time.  The  first  to  reach  maturity,  the  “zone” 
models,  describe  compartment  fires.  Each  compartment  is  divided  into  two  spatially  homogeneous  volumes, 
a hot  upper  layer  and  a cool  lower  layer.  Mass  and  energy  balances  are  enforced  for  each  layer,  with 
additional  models  describing  other  physical  processes  appended  as  differential  or  algebraic  equations  as 
appropriate.  Examples  of  such  phenomena  include  fire  plumes,  flows  through  doors,  windows  and  other 
vents,  radiative  and  convective  heat  transfer,  and  solid  fuel  pyrolysis.  An  excellent  description  of  the  physical 
and  mathematical  assumptions  behind  the  zone  modeling  concept  is  given  by  Quintiere  [2],  who  chronicles 
developments  through  1983.  Model  development  since  then  has  progressed  to  the  point  where  documented 
and  supported  software  implementing  these  models  are  widely  available  [3]. 

The  relative  physical  and  computational  simplicity  of  the  zone  models  has  led  to  their  widespread  use  in 
the  analysis  of  fire  scenarios.  So  long  as  detailed  spatial  distributions  of  physical  properties  are  not  required, 
and  the  two  layer  description  reasonably  approximates  reality,  these  models  are  quite  reliable.  However, 
by  their  very  nature,  there  is  no  way  to  systematically  improve  them.  The  rapid  growth  of  computing 
power  and  the  corresponding  maturing  of  computational  fluid  dynamics  (CFD),  has  led  to  the  development 
of  CFD  based  “field”  models  applied  to  fire  research  problems.  Virtually  all  this  work  is  based  on  the 
conceptual  framework  provided  by  the  Reynolds-averaged  form  of  the  governing  equations,  in  particular 
the  k — e turbulence  model  pioneered  by  Patankar  and  Spalding  [4],  The  use  of  CFD  models  has  allowed  the 
description  of  fires  in  complex  geometries,  and  the  incorporation  of  a wide  variety  of  physical  phenomena. 
However,  these  models  have  a fundamental  limitation  for  fire  applications  - the  averaging  procedure  at 
the  root  of  the  model  equations.  The  k — e model  was  developed  as  a time-averaged  approximation  to  the 
conservation  equations  of  fluid  dynamics.  While  the  precise  nature  of  the  averaging  time  is  not  specified,  it  is 
clearly  long  enough  to  require  the  introduction  of  large  eddy  transport  coefficients  to  describe  the  unresolved 
fluxes  of  mass,  momentum  and  energy.  This  is  the  root  cause  of  the  smoothed  appearance  of  the  results  of 
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even  the  most  highly  resolved  fire  simulations.  The  smallest  resolvable  length  scales  are  determined  by 
the  product  of  the  local  velocity  and  the  averaging  time,  rather  than  the  spatial  resolution  of  the  underlying 
computational  grid.  This  property  of  the  k — £ model  is  typically  exploited  in  numerical  computations  by 
using  implicit  numerical  techniques  to  take  large  time  steps. 

Unfortunately,  the  evolution  of  large  eddy  structures  characteristic  of  most  fire  plumes  is  lost  with 
such  an  approach,  as  is  the  prediction  of  local  transient  events.  It  is  sometimes  argued  that  the  averaging 
process  used  to  define  the  equations  is  an  “ensemble  average”  over  many  replicates  of  the  same  experiment 
or  postulated  scenario.  However,  this  is  a moot  point  in  fire  research  since  neither  experiments  nor  real 
scenarios  are  replicated  in  the  sense  required  by  that  interpretation  of  the  equations.  The  application  of 
“Large  Eddy  Simulation”  (LES)  techniques  to  fire  is  aimed  at  extracting  greater  temporal  and  spatial  fidelity 
from  simulations  of  fire  performed  on  the  more  finely  meshed  grids  allowed  by  ever  faster  computers.  The 
phrase  LES  refers  to  the  description  of  turbulent  mixing  of  the  gaseous  fuel  and  combustion  products  with 
the  local  atmosphere  surrounding  the  fire.  This  process,  which  determines  the  burning  rate  in  most  fires  and 
controls  the  spread  of  smoke  and  hot  gases,  is  extremely  difficult  to  predict  accurately.  This  is  true  not  only 
in  fire  research  but  in  almost  all  phenomena  involving  turbulent  fluid  motion.  The  basic  idea  behind  the 
LES  technique  is  that  the  eddies  that  account  for  most  of  the  mixing  are  large  enough  to  be  calculated  with 
reasonable  accuracy  from  the  equations  of  fluid  dynamics.  The  hope  (which  must  ultimately  be  justified  by 
appeal  to  experiments)  is  that  small  scale  eddy  motion  can  either  be  crudely  accounted  for  or  ignored. 

The  equations  describing  the  transport  of  mass,  momentum,  and  energy  by  the  fire  induced  flows  must 
be  simplified  so  that  they  can  be  efficiently  solved  for  the  fire  scenarios  of  interest.  The  general  equations  of 
fluid  dynamics  describe  a rich  variety  of  physical  processes,  many  of  which  have  nothing  to  do  with  fires. 
Retaining  this  generality  would  lead  to  an  enormously  complex  computational  task  that  would  shed  very 
little  additional  insight  on  fire  dynamics.  The  simplified  equations,  developed  by  Rehm  and  Baum  [5],  have 
been  widely  adopted  by  the  larger  combustion  research  community,  where  they  are  referred  to  as  the  “low 
Mach  number”  combustion  equations.  They  describe  the  low  speed  motion  of  a gas  driven  by  chemical  heat 
release  and  buoyancy  forces. 

The  low  Mach  number  equations  are  solved  numerically  by  dividing  the  physical  space  where  the  fire 
is  to  be  simulated  into  a large  number  of  rectangular  cells.  Within  each  cell  the  gas  velocity,  temperature, 
etc.,  are  assumed  to  be  uniform;  changing  only  with  time.  The  accuracy  with  which  the  fire  dynamics  can 
be  simulated  depends  on  the  number  of  cells  that  can  be  incorporated  into  the  simulation.  This  number 
is  ultimately  limited  by  the  computing  power  available.  Present  day  desktop  computers  limit  the  number 
of  such  cells  to  at  most  a few  million.  This  means  that  the  ratio  of  largest  to  smallest  eddy  length  scales 
that  can  be  resolved  by  the  computation  (the  “dynamic  range”  of  the  simulation)  is  roughly  100  ~ 200. 
Unfortunately,  the  range  of  length  scales  that  need  to  be  accounted  for  if  all  relevant  fire  processes  are  to  be 
simulated  is  roughly  10*  ~ 105  because  combustion  processes  take  place  at  length  scales  of  1 mm  or  less, 
while  the  length  scales  associated  with  building  fires  are  of  the  order  of  meters  or  tens  of  meters.  The  form 
of  the  numerical  equations  discussed  below  depends  on  which  end  of  the  spectrum  one  wants  to  capture 
directly,  and  which  end  is  to  be  ignored  or  approximated. 
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2 Hydrodynamic  Model 


An  approximate  form  of  the  Navier-Stokes  equations  appropriate  for  low  Mach  number  applications  is 
used  in  the  model.  The  approximation  involves  the  filtering  out  of  acoustic  waves  while  allowing  for  large 
variations  in  temperature  and  density  [5].  This  gives  the  equations  an  elliptic  character,  consistent  with 
low  speed,  thermal  convective  processes.  The  computation  can  either  be  treated  as  a Direct  Numerical 
Simulation  (DNS),  in  which  the  dissipative  terms  are  computed  directly,  or  as  a Large  Eddy  Simulation 
(LES),  in  which  the  large  scale  eddies  are  computed  directly  and  the  sub-grid  scale  dissipative  processes 
are  modeled.  The  choice  of  DNS  vs.  LES  depends  on  the  objective  of  the  calculation  and  the  resolution 
of  the  computational  grid.  If,  for  example,  the  problem  is  to  simulate  the  flow  of  smoke  through  a large, 
multi-room  enclosure,  it  is  not  possible  to  resolve  the  combustion  and  transport  processes  directly.  However, 
for  small-scale  combustion  experiments,  it  is  possible  to  compute  the  transport  directly  and  the  combustion 
processes  to  some  extent. 

2.1  Conservation  Equations 

First,  consider  the  conservation  equations  of  mass,  momentum  and  energy  for  a thermally-expandable, 
multi-component  mixture  of  ideal  gases  [5]: 

Conservation  of  Mass 

^ + V-pu  = 0 (1) 

Conservation  of  Species 

C(p  Y,)  + V • pr,u  = V • (pD),VY,  + Wl"  (2) 

Conservation  of  Momentum 

f+V-x  (3) 

Conservation  of  Energy 

^ (p/?)  + V • p/?u  — ——  = q + V • kWT  + V ■ ^/2/(pZ))/VT/  (4) 

Note  that  the  external  force  on  the  fluid,  represented  by  the  term  f in  Eq.  (3),  consists  of  the  drag  exerted  by 
water  droplets  emanating  from  sprinklers.  The  energy  driving  the  system  is  represented  by  the  heat  release 
rate  c/"  in  Eq.  (4).  The  term  Dp/Dt  = dp/dt  + u • Vp  is  a material  derivative.  All  other  symbols  are  listed 
in  the  Nomenclature  (Section  8). 

2.2  State,  Mass  and  Energy  Equations 

The  conservation  equations  are  supplemented  by  an  equation  of  state  relating  the  thermodynamic  quantities 
density,  pressure  and  enthalpy;  p,  p and  h.  The  pressure  is  decomposed  into  three  components 

P = Po-poogz  + p (5) 

The  first  term  on  the  right  hand  side  is  the  “background”  pressure,  the  second  is  the  hydrostatic  contribution, 
and  the  third  is  the  flow-induced  perturbation  pressure.  For  most  applications,  $ is  constant.  However,  if 


^ + (ii-V)u  ) +Vp  = pg 
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the  enclosure  is  tightly  sealed,  p)  is  allowed  to  increase  (or  decrease)  with  time  as  the  pressure  within  the 
enclosure  rises  due  to  thermal  expansion  or  falls  due  to  forced  ventilation.  Also,  if  the  height  of  the  domain 
is  on  the  order  of  a kilometer,  po  can  no  longer  be  assumed  constant  and  must  be  considered  a function  of 
the  altitude  [6], 

The  purpose  of  decomposing  the  pressure  is  that  for  low-Mach  number  flows,  it  can  be  assumed  that  the 
temperature  and  density  are  inversely  proportional,  and  thus  the  equation  of  state  can  be  approximated  [5] 

po  = p = pTOt/M  (6) 


The  pressure  p in  the  state  and  energy  equations  is  replaced  by  the  background  pressure  p to  filter  out 
sound  waves  that  travel  at  speeds  that  are  much  faster  than  typical  flow  speeds  expected  in  fire  applications. 
The  low  Mach  number  assumption  serves  two  purposes.  First,  the  filtering  of  acoustic  waves  means  that  the 
time  step  in  the  numerical  algorithm  is  bound  only  by  the  flow  speed  as  opposed  to  the  speed  of  sound,  and 
second,  the  modified  state  equation  leads  to  a reduction  in  the  number  of  dependent  variables  in  the  system 
of  equations  by  one.  The  energy  equation  (4)  is  never  explicitly  solved,  but  its  source  terms  are  included  in 
the  expression  for  the  flow  divergence,  an  important  quantity  in  the  analysis  to  follow. 

A further  assumption  about  the  thermodynamic  variables  is  that  the  constant-pressure  specific  heat  of 
the  /th  species  cpj  is  assumed  to  be  independent  of  temperature.  Under  this  assumption,  the  enthalpy  can  be 
written  as: 

h = Jih,Yi  = Tjicp,iYi  (7) 

/ / 


The  specific  heat  for  each  species  can  be  expressed  in  terms  of  the  number  of  internal  degrees  of  freedom  y 
active  in  that  molecule. 


2+vA  ^ = ( y i \ 

2 ) M,  U-l  ) M, 


(8) 


If  the  ratio  of  specific  heats  y for  each  species  is  assumed  to  be  that  of  a diatomic  molecule  (v  = 5,  y = 7/5), 
the  equation  of  state  can  be  rewritten  in  the  form1 


Po(t) 


y- 1 

Y 


p /? 


(9) 


Taking  the  material  derivative  of  Eq.  (9)  and  using  the  mass  and  energy  conservation  equations,  the  diver- 
gence of  the  velocity  field,  V • u,  can  be  written  in  terms  of  the  thermodynamic  quantities 


V ■ u = - — - ( q"  + V ■ kVT  + V • y,cp  /TpDVT/  — J 00) 

ypo  \ i Y-  1 dt  ) 

If  the  enclosure  is  tightly  sealed,  the  background  pressure  p,  can  no  longer  be  assumed  constant  due  to 
the  increase  (or  decrease)  in  mass  and  thermal  energy  within  the  enclosure.  The  evolution  equation  for  the 
pressure  is  found  by  integrating  Eq.  (10)  over  the  entire  domain  Q. 


dpo 

dt 


— — - ( / q" dV  + f kVT-dS  + yf  CpiTpDVYpdS 

V \J  q JdQ  ,Jda 


ypo 

V 


u • dS 


Isa 


(11) 


where  V is  the  volume  of  the  enclosure. 


'The  basis  of  this  approximation  is  that  nitrogen  will  be  the  dominant  species  in  most  fire  scenarios. 
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2.3  The  Momentum  Equation 

The  momentum  equation  is  simplified  by  subtracting  off  the  hydrostatic  pressure  gradient  from  the  momen- 
tum equation  (3),  and  then  dividing  by  the  density  to  obtairr 

^‘  + ^V|u|2-uxto+-V/?  = -((p-p00)g  + f+ V-x)  (12) 

at  2 p p 

To  simplify  this  equation  further,  a substitution  is  made 


V^/  « ^V|u|2  + (13) 

The  basis  for  this  approximation  is  seen  in  the  evolution  equation  for  the  circulation,  obtained  by  integrating 
Eq.  (12)  over  a closed  loop  moving  with  the  fluid  (in  the  absence  of  any  external  force) 

= /-(-vP+(p-p-)g  + v^)-^x  (14) 

at  J p 

There  are  three  sources  of  vorticity:  the  baroclinic  torque  due  to  the  non-alignment  of  the  density  and 
pressure  gradients,  buoyancy  due  to  horizontal  density  gradients,  and  viscosity.  Buoyancy  is  the  dominant 
source  of  vorticity,  and  the  approximation  above  is  equivalent  to  neglecting  the  baroclinic  torque. 

Neglecting  the  baroclinic  torque  simplifies  the  elliptic  partial  differential  equation  obtained  by  taking 
the  divergence  of  the  momentum  equation 

= - — V-F  ; F = — u x (0- - ((p  — poc)g  + f-l- V- x)  (15) 


The  linear  algebraic  system  arising  from  the  discretization  of  Eq.  (15)  has  constant  coefficients  and  can  be 
solved  to  machine  accuracy  by  a fast,  direct  ( i.e . non-iterative)  method  that  utilizes  fast  Fourier  transforms. 
No-flux  or  forced-flow  boundary  conditions  are  specified  by  asserting  that 


r du„ 

dn  ~ n dt 


(16) 


where  Fn  is  the  normal  component  of  F at  the  vent  or  solid  wall,  and  di^/dt  is  the  prescribed  rate  of  change 
in  the  normal  component  of  velocity  at  a forced  vent.  Initially,  the  velocity  is  zero  everywhere.  At  open 
external  boundaries  the  pressure-like  term  !H  is  prescribed,  depending  on  whether  the  flow  is  outgoing  or 
incoming 

X = |u|2/2  outgoing  (17) 

9~(  = 0 incoming 

The  outgoing  boundary  condition  assumes  that  the  pressure  perturbation  p is  zero  at  an  outgoing  boundary 
and  that  ‘H  is  constant  along  streamlines.  The  incoming  boundary  condition  assumes  that  is  zero  infinitely 
far  away. 

The  components  of  the  viscous  stress  tensor  are  given  by 


dw,  ditj  2 du/t 
dxj  dxj  IJ  3 dxk 


(18) 


In  the  numerical  model,  there  are  two  options  for  treating  the  dynamic  viscosity  p.  For  a Large  Eddy 
Simulation  (LES)  where  the  grid  resolution  is  not  fine  enough  to  capture  the  mixing  processes  at  all  relevant 

2Note  the  use  of  the  vector  identity  (u  • V)u  = jV|u|2  - uxo. 
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scales,  a sub-grid  scale  model  for  the  viscosity  is  applied.  Following  the  analysis  of  Smagorinsky  [7],  the 
viscosity  can  be  modeled  as 


Res  = max  (,uDNS,  P (Q  A)2  |S|) 


(19) 


where  Cs  is  an  empirical  constant,  A is  a length  on  the  order  of  the  size  of  a grid  cell,  and  |S|  is  the  magnitude 
of  the  deformation  tensor 


|A|2  = 2 


+ 2 


+ 


fdu  dv\ 

vay  + fa) 


+ 


/ dv  dw 
\ dz  dv 
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(20) 


The  thermal  conductivity  and  material  diffusivity  are  related  to  the  turbulent  viscosity  by 


k 


LES  — 


Res  Cp 

Pr 


(P£>)/,LES  = 


Res 

Sc 


(21) 


The  Prandtl  number  Pr  and  the  Schmidt  number  Sc  are  assumed  to  be  constant  for  a given  scenario. 

There  have  been  numerous  refinements  of  the  original  Smagorinsky  model  [8, 9,  10],  but  it  is  difficult  to 
assess  the  improvements  offered  by  these  newer  schemes.  There  are  two  reasons  for  this.  First,  the  structure 
of  the  fire  plume  is  so  dominated  by  the  large  scale  resolvable  eddies  that  even  a constant  eddy  viscosity 
gives  results  almost  identical  to  those  obtained  using  the  Smagorinsky  model  [11].  Second,  the  lack  of 
precision  in  most  large  scale  fire  test  data  makes  it  difficult  to  assess  the  relative  accuracy  of  each  model. 
The  Smagorinsky  model  with  constant  Q produces  satisfactory  results  for  most  large  scale  applications 
where  boundary  layers  are  not  well  resolved. 

For  a Direct  Numerical  Simulation  (DNS),  the  viscosity,  thermal  conductivity  and  material  diffusivity 
are  approximated  from  kinetic  theory.  The  viscosity  of  the  / th  species  is  given  by 


R 


26.69  x 10 -1{M,T)X2 

GjQ.v 


kg 
m s 


(22) 


where  C/  is  the  Lennard-Jones  hard-sphere  diameter  (4)  and  £2,  is  the  collision  integral,  an  empirical  func- 
tion of  the  temperature  T.  The  thermal  conductivity  of  the  /th  species  is  given  by 


, _ R cpj  W 
1 ~ Pr  m K 


(23) 


where  the  Prandtl  number  Pr  is  0.7.  The  viscosity  and  thermal  conductivity  of  a gas  mixture  are  given  by 


R ns  — ^ k/  r ; kDNS  — y,  Y/  kj 


(24) 


The  binary  diffusion  coefficient  of  the  /th  species  diffusing  into  the  mth  species  is  given  by 


Dim  — 


2.66  x 10 -7r3/2 


nr 

s 


(25) 


where  Mlm  = 2(1 /M/  + 1 alm  = (c/  + om)/2,  and  Q.D  is  the  diffusion  collision  integral,  an  empirical 

function  of  the  temperature  T [12].  It  is  assumed  that  nitrogen  is  the  dominant  species  in  any  combustion 
scenario  considered  here,  thus  the  diffusion  coefficient  in  the  species  mass  conservation  equations  is  that  of 
the  given  species  diffusing  into  nitrogen 


( P^9)/,dns  — P Dio 


where  species  0 is  nitrogen. 


(26) 
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3 Combustion 


In  the  model,  the  combustion  processes  are  modeled  in  two  ways.  In  a DNS  calculation  where  the  diffu- 
sion of  fuel  and  oxygen  can  be  modeled  directly,  a global  one-step,  finite-rate  chemical  reaction  is  usually 
used.  In  cases  where  the  computational  grid  is  not  fine  enough  to  resolve  the  diffusion  of  fuel  and  oxygen, 
Lagrangian  particles  or  “thermal  elements”  are  used  to  introduce  the  thermal  energy  of  the  fire  into  the 
calculation,  and  also  to  visualize  the  movement  of  the  smoke  and  hot  gases.  The  thermal  elements  carry 
the  heat  released  by  the  fire,  providing  a self-consistent  description  of  the  smoke  transport  at  all  resolvable 
length  and  time  scales. 

The  terms  LES  and  DNS  refer  to  solution  methodologies  of  the  equations  of  fluid  dynamics.  The  two 
relatively  simple  combustion  models  described  below  are  not  inherent  to  LES  or  DNS.  In  fact,  the  terms 
“Direct  Numerical  Simulation”  and  “Large  Eddy  Simulation”  are  more  difficult  to  define  when  used  in  the 
context  of  a reacting  flow  calculation. 


3.1  Thermal  Elements 


In  an  LES  calculation,  combustion  is  occuring  at  length  scales  well  below  the  resolution  limits  of  the  un- 
derlying numerical  grid,  thus  the  mixing  of  fuel  gases  and  air  cannot  be  calculated  directly.  Instead,  the 
fire  is  represented  by  discrete  Lagrangian  particles  (or  thermal  elements)  that  originate  at  solid  surfaces  and 
release  heat  at  a specified  rate.  Surfaces  heat  up  due  to  both  convective  and  radiative  heat  transfer  from  some 
external  source,  like  an  ignitor.  When  a surface  heats  up  to  its  prescribed  ignition  temperature,  thermal  ele- 
ments are  ejected  from  the  surface,  and  burned  at  a prescribed  rate.  The  thermal  elements  are  introduced  at 
a rate  of  n"  particles  per  unit  time  per  unit  area  with  a small  initial  velocity.  This  initial  velocity  is  a function 
of  the  mass  loss  rate  per  unit  area  of  the  fuel  bed,  which  can  be  obtained  from  the  heat  release  rate  per  unit 
area  c/L  the  density  of  the  fuel  gases  p/  , and  the  heat  of  combustion  AH 


un  = 


JL 

P/A H 


(27) 


The  heat  release  rate  of  a single  thermal  element  is  given  by 


qf  1 

^ = vTb 


(28) 


where  tp  is  the  burn-out  time  of  the  thermal  element,  and  ft  is  the  time  the  element  is  ejected  from  the  burning 
surface.  A fraction  of  the  energy  assigned  to  each  thermal  element  is  emitted  as  radiation  and  potentially  re- 
absorbed by  surrounding  smoke  and  hot  gases.  The  volumetric  heat  release  rate  term  in  the  energy  equation 
is  thus  given  by 


q = 


X/ff/yO  %r) 

6x  8y  8z 


+ K(X)X 


%r  Qp.i 


-jK(l)dl 


\ 47t|x/v--x|2 


(29) 


The  first  summation  is  over  all  active  thermal  elements  in  the  grid  cell  whose  volume  is  8xSy8z.  The 
second  summation  is  over  all  other  active  thermal  elements  outside  of  the  grid  cell  whose  radiative  energy 
is  potentially  re-absorbed  by  the  soot  contained  within  the  grid  cell.  The  fraction  of  the  elements’  energy 
converted  into  thermal  radiation  is  constant  and  given  by  yc.  The  position  of  the  ith  element  is  Xpj,  the  heat 
release  rate  is  qr and  the  segment  of  the  line  connecting  the  points  and  x is  given  by  dl.  Since  there  are 
hundreds  of  thousands  of  thermal  elements  in  a typical  calculation,  the  summation  is  made  over  a sampling 
of  the  elements  that  are  still  burning.  The  absorption  coefficient  k is  given  by 


k=  1186^  nr1 

P soot 


(30) 
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where  pT9  is  the  mass  of  soot  per  unit  volume  and  $nor  is  the  density  of  soot  itself  [13].  More  simply, 
pYs/Psoot  = fv  the  soot  volume  fraction. 

This  simple  ray-tracing  technique  to  compute  the  transport  of  radiant  energy  from  the  fire  to  its  sur- 
roundings works  well  in  cases  where  the  fire  itself  is  the  source  of  the  radiation.  However,  once  walls 
and  the  smoke  attain  temperatures  above  40(TC,  the  fire  is  no  longer  the  sole  source  of  radiation  and  the 
use  of  ray-tracing  become  more  expensive  and  cumbersome.  A better  methodology  would  be  to  solve  the 
governing  equation  for  thermal  radiation  directly.  An  effort  is  underway  to  implement  this  solver. 

The  heat  release  rate  per  unit  area  cjj  is  prescribed  by  the  user  and  is  based  on  experimental  measure- 
ments for  a configuration  resembling  as  much  as  possible  that  being  modeled.  Typically,  a calorimetry 
experiment  will  produce  a time  history  of  the  total  heat  release  rate  of  the  fuel  array.  Dividing  the  heat 
release  rate  by  the  estimated  area  of  the  burning  surface  yields  a time  history  of  the  heat  release  rate  per  unit 
area  c/'f  {t).  Often  this  function  can  be  approximated  as  a constant,  but  a time-dependent  function  can  also 
be  used  in  the  calculation. 

The  burn-out  time  tb  is  obtained  from  the  plume  correlation  of  Baum  and  McCaffrey  [14].  It  is  assumed 
that  a thermal  element  bums  out  somewhere  in  the  intermittent  region  of  the  plume,  1 .32  D <z<  3.30  D*, 
where  z is  the  height  above  the  fuel  surface,  Et  = {Q/ cppaoTaay/g)2l5  is  the  characteristic  diameter  of  the 
fire,  and  Q is  the  total  heat  release  rate  of  the  fire.  An  estimate  of  the  burn-out  time  can  be  made 


r 

.32D*  dz  ,3.300* 

— ; — - <tb< 

dz 

Jo 

w(z)  Jo 

w(z) 

where 

r 2.18^ 

(~)  { 2.45vW 

(z  < 1.32  D*) 

flame  region 

(1.32  D*  < z<  330  D* 

‘ ) intermittent  region 

(31) 


(32) 


The  bum-out  time  falls  somewhere  between  1 .0 5^/D*/g  <tb<  1 .86 \/D* / g and  is  usually  a few  tenths  of 
a second. 

The  bum-out  time  of  any  thermal  element  will  vary  based  on  the  concentration  of  oxygen  in  the  gas 
surrounding  it3.  Oxygen  is  consumed  in  any  given  grid  cell  based  on  the  amount  of  heat  generated  in  that 
grid  cell.  The  source  term  in  the  oxygen  transport  equation  becomes 


K = - 


A 


(33) 


where  A Hq2  is  the  amount  of  heat  liberated  per  unit  mass  of  oxygen  consumed  (by  default  13,100  kJ/kg). 
When  the  oxygen  mass  fraction  Yq 2 falls  to  a certain  prescribed  lower  limit,  combustion  is  assumed  to  stop, 
and  the  unbumed  fuel  associated  with  the  thermal  elements  remains  unbumed  until  more  oxygen  can  be 
found.  This  combustion  model  is  preliminary.  In  cases  where  a fire  in  a room  is  spreading  to  the  point 
of  flashover,  it  can  no  longer  be  assumed  that  there  is  a constant  flux  of  fuel  emanating  from  the  burning 
surfaces,  nor  can  it  be  assumed  that  the  fuel  even  bums.  In  simulations  approaching  flashover,  this  simple 
model  breaks  down,  and  an  improved  combustion  model  is  being  developed. 

Another  useful  input  parameter  associated  with  a given  solid  fuel  is  the  available  potential  energy  per 
unit  volume  of  the  unbumed  fuel.  If  prescribed,  the  computational  cell  representing  a piece  of  the  fuel  will 
disappear  once  the  potential  energy  has  been  liberated.  This  parameter  is  denoted  by  or  Heat  Released 
Per  Unit  Volume,  and  is  expressed  in  units  of  kJ/rrr3.  If  prescribed,  this  parameter  permits  fire  spread  through 
and  consumption  of  the  solid  fuel,  as  in  the  burning  of  a boxed  commodity  or  a wood  crib. 

?By  default,  oxygen  transport  is  not  included  in  an  LES  calculation,  and  the  burn-out  time  of  any  thermal  element  is  based  on 
standard  flame  height  correlations. 
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3.2  One-step,  Finite-rate  Combustion 


In  a DNS  calculation,  the  diffusion  of  fuel  and  oxygen  can  be  modeled  directly,  thus  it  is  possible  to  im- 
plement a relatively  simple  one-step  chemical  reaction.  Consider  the  reaction  of  oxygen  and  a hydrocarbon 
fuel 


vCxHy  C.vHv  + \’o2  O2  > Vco2  CO2  + Vh20  H2O 


(34) 


The  reaction  rate  is  given  by  the  expression 


(35) 


Suggested  values  of  /?,  E , a and  b for  various  hydrocarbon  fuels  are  given  in  Refs.  [15,  16].  It  should  be 
understood  that  the  implementation  of  any  of  these  one-step  reaction  schemes  is  still  very  much  a research 
exercise  because  it  is  not  un  ersally  accepted  that  combustion  phenomena  can  be  represented  by  such  a sim- 
ple mechanism.  Efforts  are  currently  underway  to  determine  in  what  cases  a one-step  reaction  mechanism 
provides  a valid  description  of  the  combustion  process. 

4 Thermal  Boundary  Conditions 

There  are  four  types  of  thermal  boundary  conditions  at  solid  surfaces:  adiabatic,  fixed  temperature,  thermally- 
thin  or  thermally-thick.  The  choice  of  boundary  condition  depends  on  the  surface  material  and  whether  or 
not  it  plays  a role  in  the  fire.  The  simplest  boundary  condition  is  adiabatic,  that  is,  no  heat  transfer  at  all.  In 
this  case,  the  solid  is  assumed  to  be  at  the  same  temperature  as  the  surrounding  gas.  The  next  type  of  bound- 
ary condition  is  where  a temperature  is  prescribed  at  the  solid  surface.  If  the  surface  material  is  assumed  to 
be  thermally-thick,  the  one-dimensional  heat  conduction  equation  is  applied  in  the  direction  n normal  to  the 
solid  surface 


(36) 


where  p„  cs  and  ks  are  the  density,  specific  heat  and  conductivity  of  the  material;  and  q",  c('rr  are  the 


convective,  radiative  and  re-radiative  fluxes  at  the  surface.  If  the  surface  material  is  assumed  to  be  thermally- 
thin,  its  temperature  is  governed  by  its  density,  specific  heat  and  thickness  8 


(37) 


In  this  case,  the  individual  values  of  the  parameters  p,  cs  and  8 are  not  as  important  as  their  product. 


The  heat  fluxes  to  the  surface  for  either  a thermally-thick  or  thermally-thin  material  consist  of  gains  and 
losses  from  convection  and  radiation.  The  convective  heat  flux  to  the  surface  IJ.  is  either  obtained  directly 
from  the  gradient  of  the  gas  temperature  at  the  boundary  in  a DNS  calculation 


(38) 


or  it  is  obtained  from  a correlation  of  the  form 


q"  = h AT  W/m2  ; h = C\AT\^  W/m2/K 


(39) 


in  an  LES  calculation.  Here  AT  is  the  difference  between  the  wall  and  gas  temperature,  and  C is  an  empirical 


constant  with  a default  value  of  1 .43  for  a horizontal  surface  and  0.95  for  a vertical  surface  [17], 
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The  radiative  heat  flux  to  the  surface,  4',  is  calculated  based  on  the  assumption  that  a prescribed  fraction 
of  the  heat  released  from  the  burning  thermal  elements  is  radiated  away,  and  this  energy  is  absorbed  by  the 
surrounding  surfaces,  as  well  as  by  the  smoke-laden  gas.  At  a given  point  % on  a solid  surface,  the  radiative 
flux  is  given  as 


■ II  I A A X'  tfpJ 

qr  = 2,  COS  <>/  — : ' p 

47t|x/v -x,|- 


(40) 


xp  j is  the  position  and  qp4  is  the  heat  release  rate  of  the  ith  thermal  element,  dl  is  an  element  of  the  line 
segment  connecting  the  points  and  xs , and  (}>,  is  the  angle  formed  by  the  normal  to  the  surface  and  the 
vector  xpj  — xs. 

Energy  from  the  heated  wall  is  lost  as  thermal  radiation  according  to 


•// 

q,  t 


(41) 


where  o is  the  Stefan-Boltzmann  constant  a = 5.67  x 1CT8  W/ur/K4,  and  £ is  the  emissivity  of  the  surface. 
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5 Sprinklers 


Simulating  the  effects  of  a sprinkler  spray  involves  a number  of  pieces:  predicting  activation,  computing  the 
droplet  trajectories  and  tracking  the  water  as  it  drips  onto  the  burning  commodity. 


5.1  Sprinkler  Activation 

The  temperature  of  the  sensing  element  of  a given  sprinkler  is  estimated  from  the  differential  equation  put 
forth  by  Heskestad  and  Bill  [18],  with  the  addition  of  several  terms  to  account  for  radiative  heating  and 
cooling  by  water  droplets  in  the  gas  stream  from  previously  activated  sprinklers 


dT, 

dt 


C G> 

(T,-Tm) — P|  u| 

RTIV  ’ RTF  1 


(42) 


Here  7/  is  the  link  temperature,  TK  is  the  gas  temperature  in  the  neighborhood  of  the  link,  %,  is  the  tempera- 
ture of  the  sprinkler  mount,  and  (3  is  the  volume  fraction  of  (liquid)  water  in  the  gas  stream.  The  sensitivity 
of  the  detector  is  characterized  by  the  value  of  RTI.  The  amount  of  heat  conducted  away  from  the  link  by  the 
mount  is  indicated  by  the  “C-Factor”,  C.  The  constant  G has  been  empirically  determined  by  DiMarzo  [19] 
to  be  6 x 106  K/(m/s)2,  and  its  value  is  relatively  constant  for  different  types  of  sprinklers. 


5.2  Sprinkler  Droplet  Size  Distribution 

Once  activation  is  predicted,  a sampled  set  of  spherical  water  droplets  is  tracked  from  the  sprinkler  to  either 
the  floor  or  the  burning  commodity.  In  order  to  compute  the  droplet  trajectories,  the  initial  size  and  velocity 
of  each  droplet  must  be  prescribed.  This  is  done  in  terms  of  random  distributions.  The  initial  droplet 
size  distribution  of  the  sprinkler  spray  is  expressed  in  terms  of  its  Cumulative  Volume  Fraction  (CVF), 
a function  that  relates  the  fraction  of  the  water  volume  (mass)  transported  by  droplets  less  than  a given 
diameter.  Researchers  at  Factory  Mutual  have  suggested  that  the  CVF  for  an  industrial  sprinkler  may  be 
represented  by  a combination  of  log-normal  and  Rosin-Rammler  distributions  [20] 


{d  < dm) 
(' dm  < d) 


(43) 


where  dm  is  the  median  droplet  diameter  (i.e.  half  the  mass  is  carried  by  droplets  with  diameters  of  di  or 
less),  and  y and  a are  empirical  constants  equal  to  about  2.4  and  0.6,  respectively.  The  median  drop  diameter 
is  a function  of  the  sprinkler  orifice  diameter,  operating  pressure,  and  geometry.  Research  at  Factory  Mutual 
has  yielded  a correlation  for  the  median  droplet  diameter  [21] 


d/n  „ r _ i 

— oc  We  3 

D 


(44) 


where  D is  the  orifice  diameter  of  the  sprinkler.  The  Weber  number,  the  ratio  of  inertial  forces  to  surface 
tension  forces,  is  given  by 


Ph  U2P 
Gw 


(45) 


where  pw  is  the  density  of  water,  U is  the  water  discharge  velocity,  and  q,  is  the  water  surface  tension 
(72.8  x 1(F3  N/m  at  20°C).  The  discharge  velocity  can  be  computed  from  the  mass  flow  rate,  which  is  a 
function  of  the  sprinkler’s  operating  pressure  and  K-Factor.  FM  reports  that  the  constant  of  proportionality  in 
Eq.  (44)  appears  to  be  independent  of  flow  rate  and  operating  pressure.  Three  different  sprinklers  were  tested 
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in  their  study  with  orifice  diameters  of  16.3  mm,  13.5  mm,  12.7  mm  and  the  constants  were  approximately 
4.3,  2.9,  2.3,  respectively.  The  strike  plates  of  the  two  smaller  sprinklers  were  notched,  while  that  of  the 
largest  sprinkler  was  not  [21  ]. 

In  the  numerical  algorithm,  the  size  of  the  sprinkler  droplets  are  chosen  to  mimic  the  Rosin-Rammler/log- 
normal  distribution.  A Probability  Density  Function  (PDF)  for  the  droplet  diameter  is  defined 


f(d)  = 


F'(d) 

cF 


(46) 


Droplet  diameters  are  randomly  selected  by  equating  the  Cumulative  Number  Fraction  of  the  droplet  distri- 
bution with  a uniformly  distributed  random  variable  U 


U(d)=  f(d')dd' 
Jo 


(47) 


Figure  1 displays  typical  Cumulative  Volume  Fraction  and  Cumulative  Number  Fraction  functions. 


Figure  1:  Cumulative  Volume  Fraction  and  Cumulative  Number  Fraction  functions  of  the  droplet 
size  distribution  from  a typical  industrial  scale  sprinkler.  The  median  diameter  dn  is  1 mm,  a = 0.6 

and  y - 2.43. 

Every  droplet  from  a given  sprinkler  is  not  tracked.  Instead,  a sampled  set  of  the  droplets  is  tracked. 
Typically,  1 ,000  particles  per  sprinkler  per  second  are  tracked.  The  total  number  of  droplets  represented 
by  each  computed  droplet  is  given  by  /(/?ra^),  where  mw  is  the  mass  flow  rate  of  water  from  a single 
sprinkler,  n is  the  number  of  droplets  tracked  per  sprinkler  per  second,  and is  the  average  mass  of  a 
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droplet.  The  average  mass  of  a droplet  can  be  expressed  in  terms  of  the  diameter  PDF 


md  - 


F'(d') 

~d«~ 


dd! 


(48) 


where  pu  is  the  density  of  water.  The  number  of  droplets  tracked  per  sprinkler  per  second,  n,  is  controlled 
by  the  user. 


5.3  Sprinkler  Droplet  Trajectory  in  Air 

For  a sprinkler  spray,  the  force  term  f in  Eq.  (3)  represents  the  momentum  transferred  from  the  water  droplets 
to  the  gas.  It  is  obtained  by  summing  the  force  transferred  from  each  droplet  in  a grid  cell  and  dividing  by 
the  cell  volume 

f=  1 IpQ^(ut/-u)|u,j-u| 

2 6a  8y  8z  1 ’ 

where  Cd  is  a drag  coefficient,  is  the  velocity  of  the  droplet,  u is  the  velocity  of  the  gas,  p is  the  density 
of  the  gas,  and  8.vSv8z  is  the  volume  of  the  grid  cell.  The  trajectory  of  an  individual  droplet  is  governed  by 
the  equation 

d 1 , 

— (mdud)  = mdg  - -pCdKr%{ud  - u)\ud  -u\  (50) 

where  md  is  the  mass  of  the  droplet.  The  drag  coefficient  is  a function  of  the  local  Reynolds  number 


CD 


Re 


( 24/Re  Re  < 1 

l 24(1  +0.1 5 Re0-687) /Re  1<  Re  <1000 
[ 0.44  1 000  < Re 

P |u<:/  — ul  2rd 
d 


(51) 

(52) 


where  /u  is  the  dynamic  viscosity  of  air.  The  sprinkler  spray  droplet  temperature  y and  mass  md  are  governed 
by  the  following  equations 


dTd  _ Ashd(Tg  — Td ) 
dt  cwmd 

dmd  _ Ashd(Tg  — Td) 
dt  hy 


Td  < Ty 
Td  = Ty 


(53) 

(54) 


where  cM  is  the  specific  heat  of  water,  hy  is  the  energy  of  vaporization  of  water,  A = 4tc ry  is  the  surface  area 
of  the  droplet,  Tg  is  the  gas  temperature,  Tv  is  the  vaporization  temperature  of  water,  ly  is  the  heat  transfer 
coefficient,  given  by 


hd  = 


Nuk 
2 rd 


(55) 


Nu  is  the  Nusselt  number  given  by 


Nu  — 2 + 0.6  Re2  Pr3 


(56) 


A good  approximation  for  the  Prandtl  number  is  about  0.7,  and  k is  the  thermal  conductivity  of  air. 
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5.4  Sprinkler  Droplet  Transport  on  a Surface 


When  a water  droplet  hits  a solid  horizontal  surface,  it  is  assigned  a random  horizontal  direction  and  moves 
at  a fixed  velocity  until  it  reaches  the  edge,  at  which  point  it  drops  straight  down  at  the  same  fixed  velocity. 
This  “dripping”  velocity  has  been  measured  to  be  on  the  order  of  0.5  m/s  [22],  Both  the  temperature  of  non- 
burning  surfaces  and  the  heat  release  rate  of  burning  surfaces  are  affected  by  the  droplets.  At  non-burning 
surfaces,  water  absorbs  heat  from  the  hot  surface  and  radiant  energy  from  the  fire.  The  equations  are  slightly 
different  than  those  governing  heat  transfer  to  the  airborne  drops. 


dTd  __  0.5As(hx(Ts  — Tcl)  + q") 
dt  cwmd 

dmd  _ 0.5As(hs(Ts  — Tcl)  + q") 
dt  hv 


Td  < Ty 
Td  = Ty 


(57) 

(58) 


where  Ts  is  the  temperature  of  the  surface,  cjj.  is  the  incoming  radiant  energy  flux,  and  4 is  the  heat  transfer 
coefficient  between  the  solid  surface  and  the  water  droplet,  assumed  to  be  constant  at  3000  W/m/K  [22], 


5.5  Fire  Suppression  by  Water 

The  above  two  sections  describe  heat  transfer  from  a droplet  of  water  to  a hot  gas,  a hot  solid,  or  both. 
Although  there  is  some  uncertainty  in  the  values  of  the  respective  heat  transfer  coefficients,  the  fundamental 
physics  are  fairly  well  understood.  However,  when  the  water  droplets  encounter  burning  surfaces,  simple 
heat  transfer  correlations  become  more  difficult  to  apply.  The  reason  for  this  is  that  the  water  is  not  only 
cooling  the  surface  and  the  surrounding  gas,  but  it  is  also  changing  the  pyrolysis  rate  of  the  fuel.  If  the 
surface  of  the  fuel  is  planar,  it  is  possible  to  characterize  the  decrease  in  the  pyrolysis  rate  as  a function  of 
the  decrease  in  the  total  heat  feedback  to  the  surface.  Unfortunately,  most  fuels  of  interest  in  fire  applications 
are  multi-component  solids  with  complex  geometry  at  scales  unresolvable  by  the  computational  grid. 

To  date,  most  of  the  work  in  this  area  has  been  performed  at  Factory  Mutual.  An  important  paper  on 
the  subject  is  by  Yu  et  al.  [23].  The  authors  consider  dozens  of  rack  storage  commodity  fires  of  different 
geometries  and  water  application  rates,  and  characterize  the  suppression  rates  in  terms  of  a few  global 
parameters.  Their  analysis  yields  an  expression  for  the  total  heat  release  rate  from  a rack  storage  fire  after 
sprinkler  activation 

Q = Go  e~k (?~?o)  (59) 

where  Go  is  the  total  heat  release  rate  at  the  time  of  application  &,  and  k is  a fuel-dependent  constant.  For 
the  FMRC  Standard  Plastic  commodity  k is  given  as 

k — 0.716  /<.  - 0.0 1 3 1 s“ 1 (60) 

where  m"  is  the  flow  rate  of  water  impinging  on  the  box  tops,  divided  by  the  area  of  exposed  surface  (top 
and  sides).  It  is  expressed  in  units  of  kg/rrr/s.  For  the  Class  II  commodity,  k is  given  as 

k = 0.536  m"  - 0.0040  s“‘  (61 ) 

Unfortunately,  this  analysis  is  based  on  global  water  flow  and  burning  rates.  Equation  (59)  accounts  for 
both  the  cooling  of  non-burning  surfaces  as  well  as  the  decrease  in  heat  release  rate  of  burning  surfaces.  In 
the  FDS  model,  the  cooling  of  unbumed  surfaces  and  the  reduction  in  the  heat  release  rate  are  computed 
locally,  thus  it  is  awkward  to  apply  a global  suppression  rule.  However,  the  exponential  nature  of  suppression 
by  water  is  observed  both  locally  and  globally,  thus  it  is  assumed  that  the  local  heat  release  rate  of  the  fuel 
can  be  expressed  in  the  form  [22] 

«)(»)  = «/,o(0  + (62) 
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Here  g"0(r)  is  the  heat  release  rate  per  unit  area  of  the  fuel  when  no  water  is  applied  and  /f  and  k2  are 
functions  of  the  local  water  mass  per  unit  area,  n\[„  expressed  in  units  of  kg/nr. 

k\  = a\  m"  s_l  (63) 

k2  = a2  m".  + b2  s-1  (64) 

The  linear  term  in  Eq.  (62)  is  based  on  the  observation  that  for  a boxed  commodity,  it  is  possible  for  the 
local  heat  release  rate  to  increase  as  the  fire  burns  into  the  box  and  is  protected  from  the  water  droplets 
by  material  overhead,  thus  often  a gradual  increase  in  the  heat  release  rate  is  observed  following  the  initial 
decrease  after  water  is  applied. 

To  develop  the  suppression  model  for  the  FMRC  Standard  Plastic  commodity,  19  experiments  were 
conducted  at  UL  under  a 2 MW  calorimeter  [22],  These  experiments  were  designed  as  small-scale  RDD 
(Required  Delivered  Density)  tests.  The  fuel/sprinkler  arrangement  consisted  of  four  boxes  of  the  FMRC 
Plastic  Commodity.  The  boxes  were  stacked  two  high.  The  two  stacks  were  positioned  15  cm  (6  in)  apart, 
the  same  separation  that  is  commonly  used  in  full-scale  tests.  A water  applicator  was  positioned  above  the 
boxes  to  deliver  a uniform  water  flux  onto  the  tops  of  the  boxes.  The  applicator  consisted  of  four  nozzles 
that  were  60  cm  (2  ft)  apart  and  30  cm  ( 1 ft)  above  the  plane  of  the  box  tops.  Several  nozzle  sizes  were  used, 
depending  on  the  desired  water  flow.  Table  4 lists  the  average  water  application  rate  per  unit  area  and  the 
time  of  water  application.  The  time  of  water  application  was  varied  from  30  s to  200  s.  The  water  flux  at  the 
box  top  was  varied  from  0.03  kg/rrr/s  (0.05  gpm/ft2)  to  0.66  kg/m2/s  (0.97  gpm/ft2).  The  ignition  source 


Table  1:  Time  and  Rate  of  Water  Application 


Test 

No. 

Application 
Time  (s) 

Total  Water  Flow 

Average  Water  Flux 

(L/s) 

(gpm) 

(L/nr/s) 

(gpm/ft-) 

1 

380 

0.98 

15.5 

0.66 

0.97 

2 

470 

0.57 

9.0 

0.38 

0.56 

3 

65 

0.41 

6.5 

0.28 

0.41 

4 

106 

0.41 

6.5 

0.28 

0.41 

5 

115 

0.11 

1.8 

0.074 

0.11 

6 

122 

0.11 

1.8 

0.074 

0.11 

7 

150 

0.079 

1.3 

0.053 

0.08 

8 

93 

0.11 

1.8 

0.074 

0.11 

9 

93 

0.21 

3.3 

0.14 

0.20 

10 

110 

0.21 

3.3 

0.14 

0.20 

11 

205 

0.21 

3.3 

0.14 

0.20 

12 

116 

0.16 

2.5 

0.11 

0.16 

13 

63 

0.16 

2.5 

0.11 

0.16 

14 

64 

0.28 

4.5 

0.19 

0.28 

15 

71 

0.079 

1.3 

0.053 

0.08 

16 

62 

0.047 

0.9 

0.032 

0.05 

17 

104 

0.047 

0.9 

0.032 

0.05 

18 

58 

0.079 

1.3 

0.053 

0.08 

19 

30 

0.079 

1.3 

0.053 

0.08 

was  a propane  igniter  that  consisted  of  two  parallel  12.5  mm  diameter  copper  tubes  each  30  cm  long. 

The  heat  release  rate  histories  for  the  experiments  and  the  simulations  are  given  in  Figs.  2-4.  The  decay, 
and  in  some  cases  re-growth,  of  the  fire  is  captured  reasonably  well  by  the  simulations.  A weakness  of  the 
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suppression  algorithm,  however,  is  its  reliance  on  5 empirical  coefficients  that  are  not  easily  measured.  It 
is  hoped  that  further  work  in  this  area  will  provide  more  insight  into  fire  suppression,  and  the  numerical 
algorithm  will  reflect  this  improved  understanding. 
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TOTAL  HEAT  RELEASE  RATE  (kW)  TOTAL  HEAT  RELEASE  RATE  (kW)  TOTAL  HEAT  RELEASE  RATE  (kW) 


Figure  2:  Simulated  (solid  lines)  and  experimental  (dashed  lines)  heat  release  rates  for  Tests  1,  3-7. 
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TOTAL  HEAT  RELEASE  RATE  (1000  BTU/min)  TOTAL  HEAT  RELEASE  RATE  (1000  BTU/mln)  TOTAL  HEAT  RELEASE  RATE  (1000  BTU/min) 


TOTAL  HEAT  RELEASE  RATE  (kW)  TOTAL  HEAT  RELEASE  RATE  (kW)  TOTAL  HEAT  RELEASE  RATE  (kW) 


Figure  3:  Simulated  (solid  lines)  and  experimental  (dashed  lines)  heat  release  rates  for  Tests  8-13. 
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TOTAL  HEAT  RELEASE  RATE  (1000  BTU/min)  TOTAL  HEAT  RELEASE  RATE  (1000  BTU/mln)  TOTAL  HEAT  RELEASE  RATE  (1000  BTU/mln) 


TOTAL  HEAT  RELEASE  RATE  (kW)  TOTAL  HEAT  RELEASE  RATE  (kW)  TOTAL  HEAT  RELEASE  RATE  (kW) 


TIME  (s) 


TIME  (s) 


TIME  (s) 


Figure  4:  Simulated  (solid  lines)  and  experimental  (dashed  lines)  heat  release  rates  for  Tests  14-19. 
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TOTAL  HEAT  RELEASE  RATE  (1000  BTU/mln)  TOTAL  HEAT  RELEASE  RATE  (1000  BTU/min)  TOTAL  HEAT  RELEASE  RATE  (1000  BTU/min) 


6 Numerical  Method 


This  section  presents  the  details  of  the  numerical  algorithm.  First  the  equations  that  are  being  solved  are 
presented.  Each  of  the  conservation  equations  emphasize  the  importance  of  the  velocity  divergence  and 
vorticity  fields,  as  well  as  the  close  relationship  between  the  thermally  expandable  fluid  equations  [5]  and  the 
Boussinesq  equations  for  which  the  authors  have  developed  highly  efficient  solution  procedures  [24,  25].  All 
spatial  derivatives  are  approximated  by  second  order  central  differences  and  the  flow  variables  are  updated 
in  time  using  an  explicit  second  order  predictor-corrector  scheme. 


6.1  Simplified  Equations 

Regardless  of  whether  one  is  performing  an  LES  or  a DNS  calculation,  the  overall  solution  algorithm  is  the 
same.  The  equations  derived  in  Section  2 that  are  to  be  solved  numerically  are  listed  again  here. 


Conservation  of  Mass 


dp 

dt 


+ u • Vp  = — pV  • u 


Conservation  of  Species 


Conservation  of  Momentum 


Divergence  Constraint 


+ u • VpF/  = — pF/V  • u + V • p DVT/  + Wj" 
at 


^ + u x (0+ W = - ((p-  po<,)g  + f+  V-  x) 
dt  p 


V u = 


y—  I 4"  + V . kWT  + V ■ ^cpJTpDWY,  ~ 

ypo  \ , y- 1 ^ 


(65) 


(66) 


(67) 


(68) 


Equation  of  State 

Po(l)  = pTXjiy//Ml  (69) 

/ 

Notice  that  the  source  terms  from  the  energy  conservation  equation  have  been  incorporated  into  the  diver- 
gence and  ultimately  are  involved  in  the  mass  conservation  equation.  The  temperature  is  found  from  the 
density  and  background  pressure  via  the  equation  of  state. 


6.2  Temporal  Discretization 

All  calculations  start  with  ambient  initial  conditions.  At  the  beginning  of  each  time  step,  the  quantities  (5, 
T”,  u",  and  p q are  known.  All  other  quantities  can  be  derived  from  them.  Note  that  the  superscript 
(n  + 1 )e  refers  to  an  estimate  of  the  value  of  the  quantities  at  the  (n  + 1 )st  time  step. 

1 . The  thermodynamic  quantities  p,  Yj,  and  p0  are  estimated  at  the  next  time  step  with  an  explicit  Euler 
step.  For  example,  the  density  is  estimated 

p(”+i),  = p”  _ 5f(u"  • Vp"  + p"V  • u")  (70) 

The  divergence  (V  • u)6'+1T  is  formed  from  these  estimated  thermodynamic  quantities.  The  normal 
velocity  components  at  boundaries  that  are  needed  to  form  the  divergence  are  assumed  known. 
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2.  A Poisson  equation  for  the  pressure  is  solved  with  a direct  solver 


S72!H"  = 


(V-u 


| (w+  I )('  


(V-u)" 


5 1 


- V • F" 


(71) 


Note  that  the  vector  F contains  the  convective,  diffusive  and  force  terms  of  the  momentum  equation. 
These  will  be  described  in  detail  below.  Then  the  velocity  is  estimated  at  the  next  time  step 


u(n+i)'  = u"-8f(F"  + W") 


(72) 


Note  that  the  divergence  of  the  estimated  velocity  field  is  identically  equal  to  the  estimated  divergence 
(V  • u)("+|)'  that  was  derived  from  the  estimated  thermodynamic  quantities.  The  time  step  is  checked 
at  this  point  to  ensure  that 


bt  < min 


/ bx  by  bz  \ 
\ u v ’ w ) 


(73) 


If  the  time  step  is  too  large,  it  is  reduced  so  that  it  satisfies  the  CFL  condition  and  the  procedure  starts 
from  the  beginning  of  the  time  step.  If  the  time  step  satisfies  the  stability  condition,  the  procedure 
continues. 


3.  The  thermodynamic  quantities  p,  Yj,  and  po  are  corrected  at  the  next  time  step.  For  example,  the 
density  is  corrected 


p«+i  — - (p n + p(«+!k  _ 5r(u("+1)f  • vp("+1)?  + p("+!kv  • u(,i+1) 


(74) 


The  divergence  (V  • u)*"+1)  is  derived  from  the  corrected  thermodynamic  quantities. 
4.  The  pressure  is  recomputed  using  estimated  quantities 

v^<«+,),  = 2(V.ur'-(V-»)f+'>— (V-u)»  „ r,„+n. 

bt 

The  velocity  is  then  corrected 


,«+ 1 


I u"  _|_  u(n+l)e  _ 


(75) 


(76) 


Note  again  that  the  divergence  of  the  corrected  velocity  field  is  identically  equal  to  the  corrected 
divergence. 


6.3  Spatial  Discretization 

Spatial  derivatives  in  the  governing  equations  are  written  as  second  order  accurate  finite  differences  on  a 
rectilinear  grid.  The  overall  domain  is  a rectangular  box  that  is  divided  into  rectangular  grid  cells.  Each  cell 
is  assigned  indices  i,  j and  k representing  the  position  of  the  cell  in  the  jc,  y and  z directions,  respectively. 
Scalar  quantities  are  assigned  in  the  center  of  each  grid  cell,  thus  $jk  is  the  density  at  the  nth  time  step 
in  the  center  of  the  cell  whose  indices  are  i,  j and  k.  Vector  quantities  like  velocity  are  assigned  at  cell 
faces,  thus  the  x component  of  velocity  u is  defined  at  the  faces  whose  normals  are  parallel  to  the  x-axis, 
the  y component  v is  defined  at  the  faces  whose  normals  are  parallel  to  the  y-axis,  and  the  z component  w is 
defined  at  the  faces  whose  normals  are  parallel  to  the  z-axis.  The  quantity  tfjk  is  the  x component  of  velocity 
at  the  forward  pointing  face  of  the  i jkth  cell;  -k  is  at  the  backward  pointing  face  of  the  ijkth  cell. 
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6.4  Large  Eddy  vs.  Direct  Numerical  Simulation 

The  major  difference  between  an  LES  and  a DNS  calculation  is  the  form  of  the  viscosity,  and  the  thermal 
and  material  diffusi vities.  For  a Large  Eddy  Simulation,  the  dynamic  viscosity  is  defined  at  cell  centers 

Hijk  = Pijk  (CjA)2  |S|  (77) 


where  Cs  is  an  empirical  constant,  A = (SxSvSz)1 , and 


\S\ 


+ 


( du  dv 
V dy  dx 


2 

+ 


2 

+ 


The  quantity  |S|  consists  of  second  order  spatial  differences  averaged  at  cell  centers.  The 
tivity  and  material  diffusivity  of  the  fluid  are  related  to  the  viscosity  by 


dw\2 
dv 


(78) 


thermal  conduc- 


kijk  — 


"pfiHijk 

Pr 


(P  D)ijk 


Sc 


(79) 


where  Pr  is  the  Prandtl  number  and  Sc  is  the  Schmidt  number,  both  assumed  constant.  Note  that  the  specific 
heat  cp.o  is  that  of  the  dominant  species  of  the  mixture.  Based  on  simulations  of  smoke  plumes,  (7  is  0.14, 
Pr  and  Sc  are  0.2.  There  is  no  rigorous  justification  for  these  choices. 

The  dynamic  viscosity,  thermal  conductivity  and  diffusion  coefficients  for  a DNS  calculation  are  defined 
at  cell  centers 


;4r 

II 

' — > 

kVliTjjk) 

(80) 

^ijk  — X ^l,ij 
1 

k HT,jk ) 

(81) 

D/jjk  — Di0(Tjjk) 

(82) 

where  the  values  for  each  individual  species  are  approximated  from  kinetic  theory  [12],  The  term  [fo  is 
the  binary  diffusion  coefficient  for  species  / diffusing  into  the  predominant  species  0,  usually  nitrogen.  It 
is  often  the  case  that  the  numerical  grid  is  too  coarse  to  resolve  steep  gradients  in  flow  quantities  when  the 
temperature  is  near  ambient.  However,  as  the  temperature  increases  and  the  diffusion  coefficients  increase 
in  value,  the  situation  improves.  As  a consequence,  there  is  a provision  in  the  numerical  algorithm  to  place 
a lower  bound  on  the  viscous  coefficients  to  avoid  numerical  instabilities  at  temperatures  close  to  ambient. 
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6.5  The  Mass  Transport  Equations 

Due  to  the  low  Mach  number  approximation  being  used  in  the  model,  the  mass  and  energy  equations  are 
combined  by  way  of  the  divergence.  The  divergence  of  the  flow  field  contains  much  of  the  fire-specific 
source  terms  described  above. 


6.5.1  Convective  and  Diffusive  Transport 


The  density  at  the  center  of  the  ijkth  cell  is  updated  in  time  with  the  following  predictor-corrector  scheme. 
In  the  predictor  step,  the  density  at  the  (n  + 1 )st  time  level  is  estimated  based  on  information  at  the  nth  level 


9 ijk  IT  jk 

bt 


+ (u-vp)"/A.  = -p;yv-u)” 


jk' 


‘jk 


(83) 


Following  the  prediction  of  the  velocity  and  background  pressure  at  the  (n  + l)st  time  level,  the  density  is 
corrected 


P ijk  2 \9ijk  P;/A 
Ut 


(/i+l). 


ijk 


’ >jk 


(84) 


The  species  conservation  equations  are  differenced  the  same  way 

m)iik 


(P  *1)%' 


6 1 


+ (u  ■ Vp Y,)1jk  = - (P Y,y;jk(V  • u)T,  + (V  • p DVYifiJk  i Wli 


I ijk 


jk 


ijk  ^ rr ijk 


(85) 


at  the  predictor  step,  and 


(p  *)$' 


(pYi)';,k+ipYi)) 


(n+l 

jk 


+ (u-Vp 


-(p>//);';+1,f(v-u);;+1^+(v.pDvy/);;;+1^+iT- 

(86) 

at  the  corrector  step. 

The  convective  terms  are  written  as  upwind-biased  differences  in  the  predictor  step  and  downwind- 
biased  differences  in  the  corrector  step.  In  the  expressions  to  follow,  the  symbol  ± means  + in  the  predictor 
step  and  — in  the  corrector  step.  The  opposite  is  true  for  =p. 


(u-Vp),7*  = 


2 

u“k  bx 

+ - 

2 

1 =F£V 

9i,j+l,k  - 9 ijk 

1 

1 ±£v„ 

2 

V,Jk  by 

2 V( 

1 ieM 

9ij,k+l  - 9 ijk 

- W;  ik 

i 

1 i eM. 

V 

9 ijk  ~ P>-  1 Jk 
bx 

9 ijk  ~~  9iJ—\,k 
k 8v 

9 ijk  — 9ij,k-\ 


+ 

+ 


bz 


bz 


(87) 


(u-VpT/W  = 


lT£„  (pT/),-+i Jk  - {9^i)ijk  , l±e« 


'ijk  ‘ 


8x 


+ 


' "i—  1 ,jk 


1 T£v  (9^l)i,j+\,k  ~ {9^l)ijk  ,l±e. 


' vi  jk 


8v 


i 


-v/j-i,* 


( P T/ ) ijk  (P^  l)i—],jk 

bx 

( P T/ ) / jk  (pT/)/,;-M 

by 


+ 


i 


lieM.  (pF/),M+i  - {9Yi)ijk  , 1 i e„.  (pT/)y*-(pT/)y,*-i 

'wijk gl 3 —Wijjc-I  — 


bz 


(88) 


Note  that  without  the  inclusion  of  the  e’s,  these  are  simple  central  difference  approximations.  The  e’s  are 
local  CFL  numbers,  eu  = ubt/bx , ev  = v8r/8v,  and  ew  = wbt/bz,  where  the  velocity  components  are  those 
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that  immediately  follow.  Their  role  is  to  bias  the  differencing  upwind.  Where  the  local  CFL  number  is  near 
unity,  the  difference  becomes  nearly  fully  upwinded.  Where  the  local  CFL  number  is  much  less  than  unity, 
the  differencing  is  more  centralized  [26]. 

The  divergence  in  both  the  predictor  and  corrector  step  is  discretized 
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The  thermal  and  material  diffusion  terms  are  pure  central  differences,  with  no  upwind  or  downwind  bias, 
thus  they  are  differenced  the  same  way  in  both  the  predictor  and  corrector  steps 
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The  temperature  is  extracted  from  the  density  via  the  equation  of  state 
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Because  only  species  1 through  N are  explicitly  computed,  the  summation  is  rewritten 
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In  isothermal  calculations  involving  multiple  species,  the  density  can  be  extracted  from  the  average  molec- 
ular weight 

p,,t  = ,/ - (95) 
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Again,  because  only  species  1 through  N are  explicitly  computed,  this  expression  can  be  written 
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6.5.2  Heat  Release  Rate  (LES) 


For  an  LES  calculation,  heat  is  added  to  the  flow  domain  through  the  use  of  heat  releasing  Lagrangian 
particles,  or  thermal  elements.  The  thermal  elements  are  ejected  from  burning  surfaces  with  a normal 
velocity  that  is  either  specified  by  the  user  or  automatically  determined  based  on  the  specified  heat  release 
rate  per  unit  area,  the  heat  of  combustion,  and  the  density  of  the  fuel  gases 
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The  thermal  elements  release  energy  at  a constant  rate.  A specified  fraction  of  the  energy  is  emitted  as 
thermal  radiation,  and  this  energy  can  be  re-absorbed  by  the  smoke-laden  gases  or  by  the  walls  if  desired. 
The  non-radiated  energy  from  the  thermal  elements  is  interpolated  on  the  computational  grid.  As  input,  the 
user  specifies  the  fraction  of  that  energy  lost  as  thermal  radiation,  and  the  burn-out  time  of  the  elements 
tt,.  The  heat  release  rate  per  unit  volume  of  the  ijkth  grid  cell  is  given  by  the  non-radiated  energy 
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where  c/l  is  the  heat  release  rate  per  unit  area  assigned  to  the  surface  from  which  the  /nth  element  originated, 
h"  is  the  number  of  thermal  elements  introduced  per  unit  time  per  unit  area  at  this  same  surface,  and  the 
summation  is  over  all  thermal  elements  within  the  grid  cell  whose  indices  are  i jk. 

If  desired,  the  radiated  fraction  of  the  energy  from  the  thermal  elements  can  be  re-absorbed  by  the 
smoke-laden  gases,  in  which  case  an  additional  contribution  to  the  heat  release  rate  per  unit  volume  in  a 
given  grid  cell  is  given  by 
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Note  that  here  the  summation  is  carried  out  over  all  (or  a sampling)  of  the  thermal  elements,  not  just  those 
within  the  i jkth  cell.  The  absorption  coefficient  k is  computed  at  the  center  of  the  grid  cell.  It  is  based  on 
the  mass  of  particulate  matter  and  the  temperature  within  that  cell 
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Here  fv  is  the  soot  volume  fraction,  given  by 
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where  pwot  is  the  density  of  the  soot,  and  m/v„  is  the  particulate  mass  carried  by  the  mth  thermal  element 
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Here  Xs  is  the  soot  yield  of  the  given  fuel,  and  ft  is  the  time  when  the  mth  thermal  element  was  introduced 
at  the  burning  surface. 

If  sprinklers  are  flowing,  the  water  droplets  can  cool  the  hot  gases.  The  heat  release  rate  per  unit  volume 
is  decreased 
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where  the  summation  is  carried  out  over  all  sprinkler  droplets  in  a grid  cell  of  volume  8a8v8z.  Here  T is  the 
gas  temperature,  7^  is  the  droplet  temperature.  Ad  = 4Ttr^  is  the  surface  area  of  a droplet  and  ty  — Nuk/2rd 
is  a heat  transfer  coefficient. 
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6.5.3  Heat  Release  Rate  (DNS) 

In  a DNS  calculation,  a one-step,  finite-rate  reaction  of  a hydrocarbon  fuel  is  assumed 


VCXH>  CyHv  + V0,  O2  > V(702  CO2  + Vh20  H2O 


(104) 


For  each  grid  cell,  at  the  start  of  a time  step  where  t = f and  Y//  H ijk  = YF(tn)  and  Y^  ljk  = Y0(tn),  the 
following  ODE  is  solved  numerically  with  a 2nd  order  Runge-Kutta  scheme 
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The  temperature  Ti]k  and  density  pljk  are  fixed  at  their  values  at  time  f and  the  ODE  is  iterated  from  f to 
tn+x  in  about  100  time  steps.  The  pre-exponential  factor  B , the  activation  energy  E , and  the  exponents  a and 
b are  input  parameters.  The  average  heat  release  rate  over  the  entire  time  step  is  given  by 
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where  5 1 = tn+l  —tn.  The  species  mass  fractions  are  adjusted  at  this  point  in  the  calculation  (before  the 
convection  and  diffusion  update) 
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6.5.4  Thermal  Boundary  Conditions 


Four  types  of  thermal  boundary  conditions  are  applied  at  solid  surfaces.  The  first,  and  simplest,  is  an 
adiabatic  boundary  condition  that  states  that  there  is  no  temperature  gradient  normal  to  the  surface.  It  is 
implemented  by  assigning  to  the  grid  cell  that  is  embedded  in  the  solid  (the  ghost  cell)  the  same  temperature 
as  the  first  cell  in  the  gas  (the  gas  cell). 

The  second  type  of  boundary  condition  is  where  the  solid  surface  has  a prescribed  temperature  (usually 
this  prescribed  temperature  is  a function  of  time). 

The  third  type  of  boundary  condition  assumes  the  solid  to  be  thermally-thin.  The  surface  temperature  is 
updated  in  time  according  to 
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where  Tw  is  the  wall  temperature,  6ts  is  the  time  step  used  when  updating  the  thermal  boundary  conditions 
(usually  greater  than  the  hydrodynamic  time  step  8f),  and  p,  cs,  8 are  the  input  density,  specific  heat  and 
thickness  of  the  wall.  In  a DNS  calculation  where  the  boundary  layer  is  resolved,  the  convective  flux  to  the 
wall  is  given  by 
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where  bn  is  the  size  of  a grid  cell  in  the  normal  direction  to  the  wall.  In  an  LES  calculation  where  the 
boundary  layer  is  not  resolved. 
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where  C is  an  empirical  coefficient  (0.95  for  vertical  surface;  1 .43  for  horizontal),  and  Jus  is  the  temperature 
of  the  gas  in  the  cell  bordering  the  wall.  The  radiative  flux  to  the  wall  is  given  by 


q'l  _ ^ Xr4/y»  COStym  e-JK(l)dl  U) 
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where  qpjn  is  the  heat  release  rate  of  the  mth  thermal  element  and  (J)„  is  the  angle  formed  by  the  normal  to 
the  surface  and  the  line  connecting  the  thermal  element  and  the  point  on  the  wall. 

The  fourth  type  of  thermal  boundary  condition  is  for  a thermally-thick  solid.  In  this  case,  a one  dimen- 
sional heat  transfer  calculation  is  performed  at  each  boundary  cell  designated  as  thermally-thick.  The  width 
of  the  solid  8 is  partitioned  into  N cells,  clustered  near  the  front  face.  The  cell  boundaries  are  located  at 
points  x, 

= (113) 

where  0 < i < N,  = /8c;,  8£,  = 8/2V,  and  0 < s < ] is  a measure  of  the  degree  of  clustering  of  the  cells 
at  the  front  face.  The  width  of  each  cell  is  8.y  = f 8^,  1 < i < N where  £(._i  — (/  — ^)8^.  The 
temperature  at  the  center  of  the  zth  cell  is  denoted  Tsp.  These  temperatures  are  updated  in  time  using  an 
implicit  Crank-Nicholson  scheme 
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for  1 <i<N.  The  boundary  condition  is  discretized 
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where  Ts  \ = (Ts  ] + Ts  q)/2  is  the  temperature  at  the  front  face.  Notice  that  the  radiative  emission  term  has 
been  linearized 
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The  wall  temperature  is  defined  Tw  = Ts  i = (Tsq  + Tsj) /2. 

Regardless  of  how  the  wall  temperature  is  determined,  there  are  two  ways  of  coupling  the  wall  temper- 
ature with  the  fluid  calculation.  Gas  phase  temperatures  are  defined  at  cell  centers;  the  wall  is  defined  at 
the  boundary  of  the  bordering  gas  phase  cell  and  a “ghost”  cell  inside  the  wall.  As  far  as  the  gas  phase  cal- 
culation is  concerned,  the  normal  temperature  gradient  at  the  wall  is  expressed  in  terms  of  the  temperature 
difference  between  the  “gas”  cell  and  the  “ghost”  cell.  The  wall  temperature  affects  the  gas  phase  calculation 
through  the  prescription  of  the  ghost  cell  temperature.  This  ghost  cell  temperature  has  no  physical  meaning 
on  its  own.  Only  the  difference  between  ghost  and  gas  cell  temperatures  matters,  for  this  defines  the  heat 
transfer  to  the  wall.  In  a DNS  calculation,  the  wall  temperature  is  assumed  to  be  an  average  of  the  ghost  cell 
temperature  and  the  temperature  of  the  first  cell  in  the  gas,  thus  the  ghost  cell  temperature  is  defined 


T ghost  = 2 Tw  — Tgas  (117) 

For  an  LES  calculation,  the  heat  lost  to  the  boundary  is  equated  with  an  empirical  expression 

Tg^-T^hos,  = c(7^  _ ^ | > _ Tw)  (118) 

where  8/z  is  the  distance  between  the  center  of  the  ghost  cell  and  the  center  of  the  gas  cell.  This  equation  is 
solved  for  Tghost,  so  that  when  the  conservation  equations  are  updated,  the  amount  of  heat  lost  to  the  wall  is 
equivalent  to  the  empirical  expression  on  the  right  hand  side.  Note  that  is  purely  a numerical  construct. 
It  does  not  represent  the  temperature  within  the  wall,  but  rather  establishes  a temperature  gradient  at  the  wall 
consistent  with  the  empirical  correlation. 
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6.5.5  Species  Boundary  Conditions 

At  solid  walls  there  is  no  transfer  of  mass,  thus  the  boundary  condition  for  the  / th  species  at  a wall  is  simply 


' I, ghost  — Y I, gas 


(1  19) 


where  the  subscripts  “ghost”  and  “gas”  are  the  same  as  above  since  the  mass  fraction,  like  temperature,  is 
defined  at  cell  centers.  At  forced  flow  boundaries  either  the  mass  fraction  Ym  or  the  mass  flux  rit{  of  species 
/ may  be  prescribed.  Then  the  ghost  cell  mass  fraction  can  be  derived  because,  as  with  temperature,  the 
normal  gradient  of  mass  fraction  is  needed  in  the  gas  phase  calculation.  For  cases  where  the  mass  fraction 
is  prescribed 

Yl, ghost  = -Ylw  — Yl,gas  ( 1 20) 

For  cases  where  the  mass  flux  is  prescribed,  the  following  equation  must  be  solved  iteratively 
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where  m'j  is  the  mass  flux  of  species  / per  unit  area,  n,  is  the  normal  component  of  velocity  at  the  wall 
pointing  into  the  flow  domain,  and  8/?  is  the  distance  between  the  center  of  the  ghost  cell  and  the  center  of 
the  gas  cell.  Notice  that  the  last  term  on  the  right  hand  side  is  subtracted  at  the  predictor  step  and  added  at 
the  corrector  step,  consistent  with  the  biased  upwinding  introduced  earlier. 


6.5.6  Density  Boundary  Condition 

Once  the  temperature  and  species  mass  fractions  have  been  defined  in  the  ghost  cell,  the  density  in  the  ghost 
cell  is  computed  from  the  equation  of  state 
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6.6  The  Momentum  Equation 


The  three  components  of  the  momentum  equation  are 
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The  spatial  discretization  of  the  momentum  equations  take  the  form 
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where  is  taken  at  center  of  cell  / jk , and  Fx,ijk  are  taken  at  the  side  of  the  cell  facing  in  the  forward  a 

direction,  and  Fy  ijk  at  the  side  facing  in  the  forward  y direction,  and  vt {jk  and  Fz  ijk  at  the  side  facing  in 
the  forward  z (vertical)  direction.  In  the  definitions  to  follow,  the  components  of  the  vorticity  (c$,  cov,co-)  are 
located  at  cell  edges  pointing  in  the  a,  y and  z directions,  respectively.  The  same  is  true  for  the  off-diagonal 
terms  of  the  viscous  stress  tensor:  \y  = xyz,  xxz  = xzx,  and  xxy  = Xyx.  The  diagonal  components  of  the  stress 
tensor  Xyy,  xxv,  and  x^;  the  external  force  components  {fx,fy,fz),  and  the  upwinding  bias  terms  £v,  and 
£w  are  located  at  the  respective  cell  faces. 
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The  variables  £„,  £,  and  £w  are  local  CFL  numbers  evaluated  at  the  same  locations  as  the  velocity  compo- 
nent immediately  following  them,  and  serve  to  bias  the  differencing  of  the  convective  terms  in  the  upwind 
direction.  The  subscript  i + l indicates  that  a variable  is  an  average  of  its  values  at  the  ith  and  the  (/+  l)th 
cell.  The  divergence  defined  in  Eq.  (89)  is  identically  equal  to  the  divergence  defined  by 
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The  equivalence  of  the  two  definitions  of  the  divergence  is  a result  of  the  form  of  the  discretized  equations, 
the  time-stepping  scheme,  and  the  direct  solution  of  the  Poisson  equation  for  the  pressure. 


6.6.1  Force  Terms 


The  external  force  term  components,  in  addition  to  including  the  effects  of  buoyancy,  may  also  include  the 
drag  force  from  sprinkler  droplets. 
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where  g = (gx,gy,gz)  Is  the  gravity  vector,  ij  is  the  radius  of  a droplet,  u = (ly,  v^,  wj)  the  velocity  of  a 
droplet,  Q the  drag  coefficient,  and  5a5v5c  the  volume  of  the  ijkth  cell.  The  summations  represent  all 
droplets  within  a grid  cell  centered  about  the  x,  y and  z faces  of  a grid  cell  respectively. 
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6.6.2  Time  Step 


The  time  step  is  determined  by  the  CFL  condition,  and  in  cases  of  high  viscosity,  a parabolic  stability 
criterion  typical  of  explicit  second  order  accurate  schemes 

8.v  8 y 5 z Pz/a-Sa-2  p,7A8y2  p,7A§r' 


8/  < min  , — ^ 


jk  vi  jk  wi  jk  jk  j'A  ^Fijk 
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The  estimated  velocities  and  are  tested  at  each  time  step  to  ensure  that  the  above 

condition  is  satisfied.  If  it  is  not,  then  the  time  step  is  set  to  0.8  of  its  allowed  maximum  value  and  the 
estimated  velocities  are  recomputed  (and  checked  again).  The  parabolic  stability  criterion  is  only  invoked 
for  a DNS  calculation. 


6.7  The  Pressure  Equation 


The  divergence  of  the  momentum  equation  yields  a Poisson  equation  for  the  pressure 

H+ijk  ~ Irtjk  + H-ijk  Hj+ijc-lHjk  + Hj-it  Hj,k+\  ~ 2%jk  + %j,k-\ 

Fx,ijk  Fx,i—l,jk  Fy  ijk  — Fy  i j—X  k Fz  i y A- 1 3 . 

= & — s ~sr — d,:Vu^ 
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The  lack  of  a superscript  implies  that  all  quantities  are  to  be  evaluated  at  the  same  time  level.  This  elliptic 
partial  differential  equation  is  solved  using  a direct  (non-iterative)  FFT-based  solver  that  is  part  of  a library 
of  routines  for  solving  elliptic  PDEs  called  CRAYFISHPAK  [27].  To  ensure  that  the  divergence  of  the  fluid 
is  consistent  with  the  definition  given  in  Eq.  (10),  the  time  derivative  of  the  divergence  is  defined 


at  the  predictor  step,  and  then 
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at  the  corrector  step.  The  discretization  of  the  divergence  was  given  in  Eq.  (89). 

Direct  Poisson  solvers  are  most  efficient  if  the  domain  is  a rectangular  region,  although  other  geometries 
such  as  cylinders  and  spheres  can  be  handled  almost  as  easily.  For  these  solvers,  the  no-flux  condition  (152) 
is  simple  to  prescribe  at  external  boundaries.  For  example,  at  the  floor,  z = 0,  the  Poisson  solver  is  supplied 
with  the  Neumann  boundary  condition 


Hj,  l ~ Hj, o 

8z 


— Fz,i  y,0 


(152) 


However,  many  practical  problems  involve  more  complicated  geometries.  For  building  fires,  doors  and 
windows  within  multi-room  enclosures  are  very  important  features  of  the  simulations.  These  elements  may 
be  included  in  the  overall  domain  as  masked  grid  cells,  but  the  no-flux  condition  (152)  cannot  be  directly 
prescribed  at  the  boundaries  of  these  blocked  cells.  Fortunately,  it  is  possible  to  exploit  the  relatively  small 
changes  in  the  pressure  from  one  time  step  to  the  next  to  enforce  the  no-flux  condition.  At  the  start  of  a time 
step,  the  components  of  the  convection/diffusion  term  F are  computed  at  all  cell  faces  that  do  not  correspond 
to  walls.  At  those  cell  faces  that  do  correspond  to  solid  walls,  prescribe 
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where  Fn  is  the  normal  component  of  F at  the  wall,  and  [3  is  a relaxation  factor  empirically  determined  to  be 
about  0.8  divided  by  the  time  step  6t.  The  asterisk  indicates  the  most  recent  value  of  the  pressure.  Obviously, 
the  pressure  at  this  particular  time  step  is  not  known  until  the  Poisson  equation  is  solved.  Equation  (153) 
asserts  that  following  the  solution  of  the  Poisson  equation  for  the  pressure,  the  normal  component  of  velocity 
un  will  be  driven  closer  to  zero  according  to 


diin 

dt 


-P«I7 
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This  is  approximate  because  the  true  value  of  the  velocity  time  derivative  depends  on  the  solution  of  the 
pressure  equation,  but  since  the  most  recent  estimate  of  pressure  is  used,  the  approximation  is  very  good. 
Also,  even  though  there  are  small  errors  in  normal  velocity  at  solid  surfaces,  the  divergence  of  each  blocked 
cell  remains  exactly  zero  for  the  duration  of  the  calculation.  In  other  words,  the  total  flux  into  a given 
obstruction  is  always  identically  zero,  and  the  error  in  normal  velocity  is  usually  at  least  several  orders  of 
magnitude  smaller  than  the  characteristic  flow  velocity.  When  implemented  as  part  of  a predictor-corrector 
updating  scheme,  the  no-flux  condition  at  solid  surfaces  is  maintained  remarkably  well. 

At  open  boundaries  (say  i = /),  ‘M  is  prescribed  depending  on  whether  the  flow  is  incoming  or  outgoing 
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where  / is  the  index  of  the  last  gas  phase  cell  in  the  x direction  and  ifjk  is  the  x component  of  velocity  at  the 
boundary.  The  value  of  H in  the  ghost  cell  is 
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6.8  Particle  Tracking 

Thermal  elements  are  introduced  into  the  flow  field  as  a means  of  introducing  heat  and  as  a way  to  visualize 
the  flow.  The  position  \p  of  each  thermal  element  is  governed  by  the  equations 


(157) 


The  thermal  element  positions  are  updated  according  to  the  same  predictor-corrector  scheme  that  is  applied 
to  the  other  flow  quantities.  Briefly,  the  position  Xp  of  a given  thermal  element  is  updated  according  to  the 
two  step  scheme 


x{pn+l)e  = x”  + 5fU"  (158) 

xj+1  = l-(xnp  + x(;+l)e+btu{n+l)^  (159) 

where  the  bar  over  the  velocity  vector  indicates  that  the  velocity  of  the  fluid  is  interpolated  at  the  element’s 
position. 
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7 Conclusion 


The  equations  and  numerical  algorithm  described  in  this  document  form  the  core  of  an  evolving  fire  model. 
As  research  into  specific  fire-related  phenomena  continues,  the  relevant  parts  of  the  model  can  be  improved. 
Because  the  model  was  originally  designed  to  analyze  industrial  scale  fires,  it  can  be  used  reliably  when  the 
fire  size  is  specified  and  the  building  is  relatively  large  in  relation  to  the  fire.  In  these  cases,  the  model  pre- 
dicts flow  velocities  and  temperatures  to  an  accuracy  of  10  to  20%  compared  to  experimental  measurements. 
Currently,  research  is  focussed  on  improving  both  the  gas  phase  and  solid  phase  descriptions  of  combustion 
in  the  model  so  that  simulations  involving  fire  growth  and  suppression,  especially  in  residential  sized  rooms, 
can  be  improved. 

Any  user  of  the  numerical  model  must  be  aware  of  the  assumptions  and  approximations  being  employed. 
There  are  two  issues  for  any  potential  user  to  consider  before  embarking  on  calculations.  First,  for  both  real 
and  simulated  fires,  the  growth  of  the  fire  is  very  sensitive  to  the  thermal  properties  of  the  surrounding 
materials.  Second,  even  if  all  the  material  properties  are  known,  the  physical  phenomena  of  interest  may  not 
be  simulated  due  to  limitations  in  the  model  algorithms  or  numerical  grid.  Except  for  those  few  materials 
that  have  been  studied  to  date  at  NIST,  the  user  must  supply  the  thermal  properties  of  the  materials,  and  then 
validate  the  performance  of  the  model  with  experiments  to  ensure  that  the  model  has  the  necessary  physics 
included.  Only  then  can  the  model  be  expected  to  predict  the  outcome  of  fire  scenarios  that  are  similar  to 
those  that  have  actually  been  tested. 
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8 Nomenclature 


A, 

B 

C 

CD 

Cs 

cr 

D 

D* 

dm 

E 

f 

g 

H 

h 

hi 

k 

M 

M, 

>K 

K 

Nu 

Pr 


P o 


q 

% 
4'r 
€ 
q rr 
€ 

Re 


Td 

RTI 

5 

Sc 

T 

t 

h 

u = («,  v,w) 
W!" 

We 

x = (x,y,z) 

Yi 

Y 


water  droplet  surface  area 

pre-exponential  factor  for  Arrhenius  reaction 

Sprinkler  C-Factor 

drag  coefficient 

Smagorinsky  constant  (LES) 

constant  pressure  specific  heat 

diffusion  coefficient 

characteristic  fire  diameter 

median  volumetric  droplet  diameter 

activation  energy 

external  force  vector  (excluding  gravity) 
acceleration  of  gravity 
total  pressure  divided  by  the  density 
enthalpy:  heat  transfer  coefficient 
enthalpy  of  /th  species 

thermal  conductivity;  suppression  decay  factor 

molecular  weight  of  the  gas  mixture 

molecular  weight  of  /th  gas  species 

water  flux  per  unit  area 

water  mass  per  unit  area 

Nusselt  number 

Prandtl  number 

pressure 

background  pressure 

pressure  perturbation 

heat  release  rate  per  unit  volume 

fire  heat  release  rate  per  unit  area 

radiative  flux  to  a solid  surface 

convective  flux  to  a solid  surface 

radiative  loss  from  a solid  surface 

water  cooling  per  unit  area 

universal  gas  constant 

Reynolds  number 

water  droplet  radius 

Response  Time  Index  of  sprinkler 

deformation  tensor 

Schmidt  number 

temperature 

time 

thermal  element  burn-out  time  (LES) 
velocity  vector 

production  rate  of  /th  species  per  unit  volume 

Weber  number 

position  vector 

mass  fraction  of  /th  species 

ratio  of  specific  heats 
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Q N N 


A H 

5 

e 

K 

V 

V, 

P 

X 


CO  = (coA-,cov,coz) 


heat  of  combustion 
wall  thickness 
emissivity 

absoiption  coefficient 
dynamic  viscosity 
stoichiometric  coefficient 
density 

viscous  stress  tensor 
radiative  loss  fraction 
smoke  or  soot  yield 
Stefan-Boltzmann  constant 
vorticity  vector 
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